{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "np.set_printoptions(precision=14, linewidth=120)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def integrate(fn, a, b, steps=5, debug=False, exact=None):\n",
    "    table = np.zeros((steps, steps), dtype=np.float64)\n",
    "    pow_4 = 4 ** np.arange(steps, dtype=np.float64) - 1\n",
    "\n",
    "    # trapezoidal rule\n",
    "    h = (b - a)\n",
    "    table[0, 0] = h * (fn(a) + fn(b)) / 2\n",
    "\n",
    "    for j in range(1, steps):\n",
    "        h /= 2\n",
    "\n",
    "        # extended trapezoidal rule\n",
    "        table[j, 0] = table[j - 1, 0] / 2\n",
    "        table[j, 0] += h * np.sum(fn(a + i * h) for i in range(1, 2 ** j + 1, 2))\n",
    "\n",
    "        # richardson extrapolation\n",
    "        for k in range(1, j + 1):\n",
    "            table[j, k] = table[j, k - 1] + (table[j, k - 1] - table[j - 1, k - 1]) / pow_4[k]\n",
    "\n",
    "    # debug\n",
    "    if debug:\n",
    "        print(table, file=sys.stderr)\n",
    "        if exact is not None:\n",
    "            errors = ['%.2e' % i for i in np.abs(table.diagonal() - exact)]\n",
    "            print('abs. error:', errors, file=sys.stderr)\n",
    "\n",
    "    return table[-1, -1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[[ 0.68393972058572  0.                0.                0.                0.              ]\n",
      " [ 0.73137025182856  0.74718042890951  0.                0.                0.              ]\n",
      " [ 0.74298409780038  0.74685537979099  0.74683370984975  0.                0.              ]\n",
      " [ 0.7458656148457   0.74682612052747  0.7468241699099   0.74682401848228  0.              ]\n",
      " [ 0.74658459678822  0.74682425743573  0.74682413322961  0.74682413264739  0.74682413309509]]\n",
      "abs. error: ['6.29e-02', '3.56e-04', '9.58e-06', '1.14e-07', '2.83e-10']\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.74682413309509432"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# integral[0, 1] of e^(-x^2)\n",
    "integrate(lambda x: np.exp(-x * x), 0, 1, debug=True, exact=0.746824132812427)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[[ 0.75              0.                0.                0.                0.              ]\n",
      " [ 0.70833333333333  0.69444444444444  0.                0.                0.              ]\n",
      " [ 0.69702380952381  0.69325396825397  0.6931746031746   0.                0.              ]\n",
      " [ 0.69412185037185  0.69315453065453  0.69314790148123  0.69314747764483  0.              ]\n",
      " [ 0.69339120220753  0.69314765281942  0.69314719429708  0.69314718307193  0.69314718191674]]\n",
      "abs. error: ['5.69e-02', '1.30e-03', '2.74e-05', '2.97e-07', '1.36e-09']\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.69314718191674496"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# ln(2)\n",
    "integrate(1..__truediv__, 1, 2, debug=True, exact=np.log(2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[[ 0.5           0.            0.            0.            0.          ]\n",
      " [ 0.3125        0.25          0.            0.            0.          ]\n",
      " [ 0.265625      0.25          0.25          0.            0.          ]\n",
      " [ 0.25390625    0.25          0.25          0.25          0.          ]\n",
      " [ 0.2509765625  0.25          0.25          0.25          0.25        ]]\n",
      "abs. error: ['2.50e-01', '0.00e+00', '0.00e+00', '0.00e+00', '0.00e+00']\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.25"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# integral[0, 1] of x^3\n",
    "integrate(lambda x: x**3, 0, 1, debug=True, exact=.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## logarithmus naturalis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ln(x):\n",
    "    if x <= 0:\n",
    "        raise ValueError()\n",
    "    m, e = np.frexp(x)\n",
    "    return integrate(1..__truediv__, 1, m) + e * 0.6931471805599453"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.99999999999789502"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ln(np.e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.69314717920314561"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ln(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1447298858493882"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ln(np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## normal distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def norm_pdf(x, mean, sd):\n",
    "    x0 = (x - mean) ** 2\n",
    "    v2 = 2 * sd ** 2\n",
    "    return np.exp(-x0 / v2) / np.sqrt(np.pi * v2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def norm_cdf(x, mean, sd):\n",
    "    return integrate(lambda x: norm_pdf(x, mean, sd), mean, x) + .5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.97500211189427477"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm_cdf(1.96, mean=0, sd=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.68268949213754959"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm_cdf(1.2, mean=1, sd=.2) - norm_cdf(0.8, mean=1, sd=.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
